import time
import matplotlib.pyplot as plt
import numpy as np
from xray_scan import XRayScan
from xray_scan_collection import XRayScanCollection
from xray_scan_collection_fft import XRayScanCollectionFFT
from xray_rec_interpolation_and_rectangle import XRayReconstruct

# ------------------各种初始参量------------------------------#
# 设定矩形区域的长和宽
rect_width = 50  # 宽度
rect_height = 50  # 高度

# 生成两个点的坐标（在矩形内部）
P1 = (10, 0)  # 第一个点的坐标 (x1, y1)
P2 = (0, 0)  # 第二个点的坐标 (x2, y2)
# 假设P1和P2的权重
P1_weight = 50
P2_weight = 30
weight_bias = 0

# 设置屏幕分辨率：将屏幕线段分为 N 等份

N = XRayReconstruct().get_all_num()  # 分辨率，即传感器数量 请始终保持为奇数 这样比较好计算
# 设置角度分辨率 也就是把180度分成多少分
S = 24

point_radius = 0

# 开关变量：是否显示分段和编号信息
show_segments = False




# 设置紫色为X射线的颜色
xray_color = 'purple'

# -----------------初始参量计算----------------------------- #

# 计算外接圆的半径（矩形对角线的一半）
circle_radius = np.sqrt((rect_width / 2) ** 2 + (rect_height / 2) ** 2)

# 计算正方形的对角线长度（与切线段等长）
line_length = np.sqrt(rect_width ** 2 + rect_height ** 2)

#XRayScan 列表
XRSCollection = XRayScanCollection()

# ------------------------------------------------------- #

# 计算与X射线发射端和屏幕相交的线段的端点
def find_intersection(line_start, line_end, slope, y_intercept):
    x1, y1 = line_start
    x2, y2 = line_end
    if x1 == x2:  # 垂直线
        x = x1
        y = slope * x + y_intercept
    else:
        intercept_line_slope = (y2 - y1) / (x2 - x1)
        intercept_line_intercept = y1 - intercept_line_slope * x1

        if intercept_line_slope == slope:  # 平行线
            return None

        x = (intercept_line_intercept - y_intercept) / (slope - intercept_line_slope)
        y = slope * x + y_intercept

    return np.array([x, y])


def generate_ray_line(ray_point, tangent_point, slope):
    y_intercept = ray_point[1] - slope * ray_point[0]
    intersect_start = find_intersection(line_start, line_end, slope, y_intercept)
    intersect_end = find_intersection(screen_line_start, screen_line_end, slope, y_intercept)
    return intersect_start, intersect_end

def find_segment(ray_end, screen_points):
    for i in range(len(screen_points) - 1):
        segment_start = screen_points[i]
        segment_end = screen_points[i + 1]
        if (segment_start[0] <= ray_end[0] <= segment_end[0]) or (segment_start[0] >= ray_end[0] >= segment_end[0]):
            return i + 1
    return None

# ----------------------------开始循环---------------------------#




for i in range(S):

    # -----------------绘制初始图像-----------------------------


    # 指定切点与x轴正半轴的夹角（角度制）
    angle_deg = 180 / S * i + 1 # 可以修改角度值 加1 为了避免0度和180度时出现问题

    angle_rad = np.deg2rad(angle_deg)  # 转换为弧度制

    print(angle_deg)

    # 计算“X射线发射端”的切点的坐标
    tangent_point_x = circle_radius * np.cos(angle_rad)
    tangent_point_y = circle_radius * np.sin(angle_rad)
    tangent_point = np.array([tangent_point_x, tangent_point_y])

    # 计算切线的方向向量（与圆的法线垂直）
    tangent_vector = np.array([-np.sin(angle_rad), np.cos(angle_rad)])
    tangent_vector /= np.linalg.norm(tangent_vector)  # 归一化

    # 计算“X射线发射端”线段的两个端点
    half_length = line_length / 2
    line_start = tangent_point - half_length * tangent_vector
    line_end = tangent_point + half_length * tangent_vector

    # 计算“屏幕”的切点（与“X射线发射端”中心对称）
    screen_tangent_point = -tangent_point  # 对称点

    # 计算“屏幕”线段的两个端点
    screen_line_start = screen_tangent_point - half_length * tangent_vector
    screen_line_end = screen_tangent_point + half_length * tangent_vector

    interval_length = line_length / N  # 每个区间的长度

    # 计算每个分割点的位置
    screen_points = [screen_line_start + i * interval_length * tangent_vector for i in range(N + 1)]

    # 计算与角度对应的斜率
    slope = np.tan(angle_rad)

    ray1_start, ray1_end = generate_ray_line(P1, tangent_point, slope)
    ray2_start, ray2_end = generate_ray_line(P2, tangent_point, slope)

    # 打印X射线在“屏幕”线段上的坐标位置
    print(f"Ray 1 intersects screen at: {ray1_end}")
    print(f"Ray 2 intersects screen at: {ray2_end}")


    # 确定每条X射线落在哪个区间



    ray1_segment = find_segment(ray1_end, screen_points)
    ray2_segment = find_segment(ray2_end, screen_points)

    print(f"Ray 1 falls into segment: {ray1_segment}")
    print(f"Ray 2 falls into segment: {ray2_segment}")



    # 初始化传感器图像
    sensor_image = np.zeros(N)

    for i in range(len(sensor_image)):
        sensor_image[i] = weight_bias

    # 找到X射线落在传感器的哪个格点上，并赋予权重
    if ray1_segment is not None:
        for i in range(-point_radius,point_radius + 1,1):
            sensor_image[ray1_segment - 1 + i] += P1_weight

    if ray2_segment is not None:
        for i in range(-point_radius, point_radius + 1, 1):
            sensor_image[ray2_segment - 1 + i] += P2_weight

    # 打印传感器图像
    # print("Sensor Image:", sensor_image)

    # 绘制传感器图像为折线图
    # plt.figure(figsize=(10, 4))
    # plt.plot(range(1, N + 1), sensor_image, 'o-', alpha=0.6)
    # plt.xlabel('Sensor Position')
    # plt.ylabel('Intensity')
    # plt.title('Sensor Image')
    # plt.grid(True, linestyle='--', alpha=0.6)
    # plt.xticks(range(1, N + 1))
    # plt.show()

    # v8 更改部分：对传感器图像做傅里叶变换并绘制功率谱
    # sensor_image_fft = np.fft.fft(sensor_image)
    # frequency = np.fft.fftfreq(N)
    # power_spectrum = np.abs(sensor_image_fft) ** 2
    #
    # positive_freq_indices = frequency >= 0
    # positive_frequencies = frequency[positive_freq_indices]
    # positive_power_spectrum = power_spectrum[positive_freq_indices]

    # plt.figure(figsize=(10, 4))
    # plt.plot(positive_frequencies, positive_power_spectrum, 'o-', alpha=0.6)
    # plt.xlabel('Frequency')
    # plt.ylabel('Power')
    # plt.title('Power Spectrum')
    # plt.grid(True, linestyle='--', alpha=0.6)
    # plt.show()

    # 对傅里叶变换结果进行反变换
    # sensor_image_ifft = np.fft.ifft(sensor_image_fft)

    # 绘制反变换后的图像
    # plt.figure(figsize=(10, 4))
    # plt.plot(range(1, N + 1), sensor_image_ifft.real, 'o-', alpha=0.6)
    # plt.xlabel('Sensor Position')
    # plt.ylabel('Intensity')
    # plt.title('Inverse FFT of Sensor Image')
    # plt.grid(True, linestyle='--', alpha=0.6)
    # plt.xticks(range(1, N + 1))
    # plt.show()

    XRS = XRayScan(angle_deg,sensor_image)
    XRSCollection.add_scan(XRS)



    # 绘制屏幕上的分割线或标记点（根据开关变量）
    if show_segments:
        # 绘制矩形区域、外接圆、两个点和切线段
        plt.figure(figsize=(10, 10))  # 放大图像
        plt.plot(P1[0], P1[1], 'ro', label='Point P1')  # 绘制第一个点
        plt.plot(P2[0], P2[1], 'bo', label='Point P2')  # 绘制第二个点

        # 绘制矩形边框
        rect = plt.Rectangle((-rect_width / 2, -rect_height / 2), rect_width, rect_height,
                             edgecolor='black', facecolor='none', linestyle='--', label='Rectangle')
        plt.gca().add_patch(rect)

        # 绘制外接圆
        circle = plt.Circle((0, 0), circle_radius, color='green', fill=False, linestyle='--',
                            label='Circumscribed Circle')
        plt.gca().add_patch(circle)
        # 绘制“X射线发射端”线段
        plt.plot([line_start[0], line_end[0]], [line_start[1], line_end[1]],
                 'm-', linewidth=2, label='X-ray Source')

        # 绘制“屏幕”线段
        plt.plot([screen_line_start[0], screen_line_end[0]], [screen_line_start[1], screen_line_end[1]],
                 'c-', linewidth=2, label='Screen')
        for i, point in enumerate(screen_points):
            plt.plot(point[0], point[1], 'k|', markersize=10, markeredgewidth=2)  # 绘制分割标记
            # 为每个传感器标记横坐标
            if i > 0:  # 跳过原点
                plt.text(point[0] + 0.2, point[1] + 0.2, f'{i}', fontsize=10, color='blue')

        # 绘制屏幕分辨率的原点
        plt.plot(screen_line_start[0], screen_line_start[1], 'go', markersize=8, label='Screen Origin')

        # 绘制与角度对应的线段
        if ray1_start is not None and ray1_end is not None:
            plt.plot([ray1_start[0], ray1_end[0]], [ray1_start[1], ray1_end[1]], color=xray_color, linestyle='--',
                     label='X-ray Ray')
        if ray2_start is not None and ray2_end is not None:
            plt.plot([ray2_start[0], ray2_end[0]], [ray2_start[1], ray2_end[1]], color=xray_color, linestyle='--')

        # 设置图形显示的范围，扩大视野
        padding = 2  # 额外的边界范围
        plt.xlim(-circle_radius - line_length / 2 - padding, circle_radius + line_length / 2 + padding)
        plt.ylim(-circle_radius - line_length / 2 - padding, circle_radius + line_length / 2 + padding)

        # 添加坐标轴
        plt.axhline(0, color='black', linewidth=0.5)
        plt.axvline(0, color='black', linewidth=0.5)

        # 添加标签和标题
        plt.xlabel('X-axis')
        plt.ylabel('Y-axis')
        plt.title('Rectangle Region with X-ray Source, Screen, and Screen Resolution Origin')
        plt.legend()
        plt.grid(True, linestyle='--', alpha=0.6)
        plt.gca().set_aspect('equal', adjustable='box')  # 保持坐标比例一致
        plt.show()

    # 显示图形

        time.sleep(1)

#画出sino图像
XRSCollection.plot_sinogram()
print(XRSCollection)
# 创建XRayScanCollectionFFT对象，并计算傅里叶变换

# 首先传入FFT类
fft_collection = XRayScanCollectionFFT(XRSCollection)

# 绘制原始数据的直角坐标图  此步骤一定要在计算FFT之前做 一定
# 这一步仅仅是绘制了图像
fft_collection.plot_original_cartesian()

# 随后计算FFT 计算过后 应该就会进行拼接
fft_collection.compute_fft()

# 绘制位于指定长方形内的点和长方形外框
aspect_ratio = 1.0  # 1:1 长宽比
fft_collection.plot_filtered_cartesian(aspect_ratio)

rectangle_data = fft_collection.get_rectangle_data(aspect_ratio)

reconstructor = XRayReconstruct()
interpolated_values_rbf = reconstructor.RBF_interpolation(rectangle_data)


reconstructor.RBF_plot(interpolated_values_rbf)



###     对 interpolated_values_rbf 进行低通滤波   ###

def lowpass_filter(interpolated_values_rbf, cutoff_frequency=0.5, filter_type='ideal'):
    """
    对插值结果进行二维低通滤波，使用方形滤波区域
    参数:
    - interpolated_values_rbf: 2D复数数组，插值后的频域数据
    - cutoff_frequency: 截止频率（0到1之间），表示保留中心区域的比例
    - filter_type: 'ideal' 或 'gaussian'，选择滤波器类型
    """
    rows, cols = interpolated_values_rbf.shape
    crow, ccol = rows // 2, cols // 2  # 中心点位置

    # 创建网格
    x = np.linspace(-1, 1, cols)
    y = np.linspace(-1, 1, rows)
    X, Y = np.meshgrid(x, y)

    # 创建滤波器
    if filter_type == 'ideal':
        # 理想低通滤波器（方形）
        H = np.zeros((rows, cols))
        mask = (np.abs(X) <= cutoff_frequency) & (np.abs(Y) <= cutoff_frequency)
        H[mask] = 1
    else:
        # 高斯低通滤波器（方形）
        sigma = cutoff_frequency / 2
        H = np.exp(-(X**2)/(2*sigma**2)) * np.exp(-(Y**2)/(2*sigma**2))

    # 应用滤波器
    filtered_values = interpolated_values_rbf * H

    return filtered_values

# 高通滤波
def high_pass_filter(data, cutoff_frequency=0.5, order=2):
    """
    应用高通滤波器到二维傅里叶空间数据

    参数:
    data: 2D numpy array, 输入的傅里叶空间数据
    cutoff_frequency: float, 截止频率 (0到1之间)
    order: int, 滤波器阶数

    返回:
    filtered_data: 经过滤波的数据
    """
    rows, cols = data.shape
    center_row, center_col = rows // 2, cols // 2

    # 创建频率网格
    y, x = np.ogrid[-center_row:rows - center_row, -center_col:cols - center_col]
    radius = np.sqrt((x / float(center_col)) ** 2 + (y / float(center_row)) ** 2)

    # 构建高通滤波器
    high_pass = 1 - 1 / (1 + (radius / cutoff_frequency) ** (2 * order))

    # 应用滤波器
    filtered_data = data * high_pass

    return filtered_data


# 使用示例：
# 理想低通滤波
filtered_ideal = lowpass_filter(interpolated_values_rbf=interpolated_values_rbf, cutoff_frequency=0.3, filter_type='ideal')

# 可视化结果
def plot_results(filtered_gaussian):
    plt.figure(figsize=(5, 5))

    # 高斯低通滤波结果
    plt.imshow(np.abs(filtered_gaussian), cmap='viridis')
    plt.title('Lowpass Filter')
    plt.colorbar()

    plt.tight_layout()
    plt.show()

filtered_ideal = np.fft.fftshift(filtered_ideal)

plot_results(interpolated_values_rbf)

interpolated_values_rbf = np.fft.ifftshift(interpolated_values_rbf)

# 进行逆变换
final_image = np.fft.ifft2(interpolated_values_rbf)

# 最终移位
final_image = np.fft.fftshift(final_image)

#调整位置
final_image = np.rot90(final_image, k=-1)
final_image = np.flip(final_image, axis=0)

# 显示结果
plt.figure(figsize=(5, 5))
plt.imshow(final_image.real)
plt.title('Properly Aligned Reconstruction')
plt.colorbar()
plt.show()


